﻿// 平差结果输出，格式ou2,ouf...
using System;
using System.Collections.Generic;
using System.Text;
using MathNet.Numerics.LinearAlgebra;
using System.IO;
using static AdjPlane.Angle_t;
using static AdjPlane.AdjCmn;

namespace AdjPlane
{
    using mat = Matrix<double>;
    using pindex = Dictionary<string, int>;

    //输出文件
    public partial class Plane
    {
        readonly int sn = 13;//点名的英文长度

        //判断是否一三象限
        protected static bool Location13(double T)
        {
            if ((T >= 0 && T <= 90) || (T >= 180 && T <= 270))
                return true;
            else
                return false;
        }

        //陀螺边相关信息
        void OutputGyro(string path)
        {
            List<string> str = new List<string>();
            int i;
            str.Add(string.Format("{0,30}", "陀螺边改正"));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0,15}{1,15}{2,15}{3,10}{4,15}", "From", "To",
                 "Gyro_Obs/dms", "V/\"", "Result/dms"));
            for(i=0;i<Gyro.Count;i++)
            {
                int ids = FindStaName(Gyro[i].beg);
                int ido = FindObsName(ids, Gyro[i].end);
                double vt = AllObs[ids].data[ido].vt;
                double gyro= D2DMS(DMS2D(this.Gyro[i].gyazi)+ vt / 3600);
                this.Gyro[i].gypost = gyro;
                str.Add(string.Format("{0,15}{1,15}{2,15}{3,10}{4,15}", this.Gyro[i].beg, this.Gyro[i].end,
                    this.Gyro[i].gyazi.ToString("f5"), vt.ToString("f2"), gyro.ToString("f5")));
            }
            str.Add("\n");
            str.Add(string.Format("{0,30}", "先验陀螺边-坐标方位角-比较"));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0,15}{1,15}{2,15}{3,15}{4,10}", "From", "To",
                   "Gyro/dms", "Azi/dms", "Dif/\""));
            foreach (GyroAzi t in Gyro)
            {
                double dx = Result[t.end].X - Result[t.beg].X;
                double dy = Result[t.end].Y - Result[t.beg].Y;
                double azi = RD2D(GetAzi(dx, dy));
                double azi_dms = D2DMS(azi);
                double da = (DMS2D(t.gyazi) - azi)*3600;
                str.Add(string.Format("{0,15}{1,15}{2,15}{3,15}{4,10}", t.beg, t.end, 
                    t.gyazi.ToString("f5"), azi_dms.ToString("f5"), da.ToString("f1")));
            }
            str.Add("\n");
            str.Add(string.Format("{0,30}", "验后陀螺边-坐标方位角-比较"));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0,15}{1,15}{2,15}{3,15}{4,10}", "From", "To",
                    "Gyro_post/dms", "Azi/dms", "Dif/\""));
            foreach (GyroAzi gy in Gyro)
            {
                double dx = Result[gy.end].X - Result[gy.beg].X;
                double dy = Result[gy.end].Y - Result[gy.beg].Y;
                double azi = RD2D(GetAzi(dx, dy));
                double azi_dms = D2DMS(azi);
                double da = (DMS2D(gy.gypost) - azi) * 3600;
                str.Add(string.Format("{0,15}{1,15}{2,15}{3,15}{4,10}", gy.beg, gy.end,
                    gy.gypost.ToString("f5"), azi_dms.ToString("f5"), da.ToString("f1")));
            }
            str.Add("\n");
            File.AppendAllLines(path, str);
        }

        /// <summary>
        /// 获得输出文件的文件名与扩展名及路径
        /// </summary>
        /// <param name="path">输入的文件</param>
        /// <param name="popt">平差参数</param>
        /// <param name="outpath">输出路径</param>
        protected void Name(Proc_t popt, ref File_t fopt)
        {
            string dir = Path.GetDirectoryName(fopt.infile);
            string p = Path.GetFileNameWithoutExtension(fopt.infile);
            string name ="",map="";
            switch (popt.WgtDef)
            {
                case Weight.Tradition: name = p + ".ou2"; break;
                case Weight.hmt1: 
                case Weight.hmt2: 
                case Weight.hmt3:
                case Weight.Helmert:
                case Weight.MINQUE:
                case Weight.LSVCE:
                case Weight.MLE:   name = p + ".ouf"; break;
                case Weight.Robust:  name = p + "_Robust.ou2"; break;
                case Weight.RobustHelmert:  name = p + "_RH.ouf"; break;
                case Weight.DataSnooping:name = p + "_DSP.ou2";break;
                default:name = p + ".ou2";break;
            }
            map = p + "2.gra";
            if (popt.WgtDef >= Weight.Helmert && popt.WgtDef <= Weight.MLE)
                map = p + "f.gra";
            fopt.oufile = dir +"\\"+ name;
            fopt.mapfile = dir +"\\"+ map;
        }
        //输出标题 
        protected void OutputTitle(Proc_t popt,string outpath)
        {
            List<string> str = new List<string>();
            switch (popt.WgtDef)
            {
                case Weight.Tradition:
                    if (popt.GyOpt.UseGyro == 1) str.Add(string.Format("{0,30}", "平面网平差结果（陀螺边约束）"));
                    else  str.Add(string.Format("{0,30}", "平面网平差结果（常规定权）")); 
                    break;
                case Weight.hmt1:
                case Weight.hmt2:
                case Weight.hmt3:
                case Weight.Helmert:  str.Add(string.Format("{0,30}", "平面网平差结果（Helmert方差分量估计）")); break;
                case Weight.MINQUE:  str.Add(string.Format("{0,30}", "平面网平差结果（MINQUE定权）")); break;
                case Weight.LSVCE:   str.Add(string.Format("{0,30}", "平面网平差结果（最小二乘估计定权）")); break;
                case Weight.MLE:  str.Add(string.Format("{0,30}", "平面网平差结果（极大似然估计定权）")); break;
                case Weight.Robust:  str.Add(string.Format("{0,30}", "平面网平差结果（抗差估计）")); break;
                case Weight.RobustHelmert:  str.Add(string.Format("{0,30}", "平面网平差结果（抗差HVCE）")); break;
            }
            str.Add("-----------------------------------------------------------------------------");
            File.WriteAllLines(outpath, str);
        }
        //输出总体信息
        protected void OutputProfile(Proc_t popt,string outpath)
        {
            int num = 30;
            List<string> str = new List<string>
            {
                string.Format("{0,30}", "平面控制网总体信息"),
                "----------------------------------------------------------------------------------",
                string.Format("{0}{1,-10}",PadLeftEx("已知点数:",num), knownPoint.Count),
                string.Format("{0}{1,-10}",PadLeftEx("未知点数:",num),ukpIndex.Count),
                string.Format("{0}{1,-10}",PadLeftEx( "测站数:",num), CZ.Count),
                string.Format("{0}{1,-10}",PadLeftEx( "方向观测数:",num),nl),
                string.Format("{0}{1,-10}",PadLeftEx( "距离观测数:",num), ns),
                string.Format("{0}{1,-10}",PadLeftEx("陀螺边数:",num),nt),
                string.Format("{0}{1,-10}",PadLeftEx("总观测数:",num),n),
                string.Format("{0}{1,-10}",PadLeftEx("未知量个数:",num),t),
                string.Format("{0}{1,-10}",PadLeftEx( "PVV:",num), PVV.ToString("f2")),
                string.Format("{0}{1,-10}",PadLeftEx( "自由度:",num), r),
                string.Format("{0}{1,-10}",PadLeftEx( "验前单位权中误差:",num), M0.ToString("f2")),
                string.Format("{0}{1,-10}",PadLeftEx( "验后单位权中误差:",num), sigma.ToString("f2")),
                "----------------------------------------------------------------------------------",
                string.Format("{0,20}{1,15}{2,10}{3,10}{4,10}", "TYPE", "PVV", "R", "D", "M")
            };
            double mL2 = PVV1 / Lri_sum;
            double mS2 = PVV2 / Sri_sum;
            str.Add(string.Format("{0,20}{1,15}{2,10}{3,10}{4,10}", "L", PVV1.ToString("f3"),
                Lri_sum.ToString("f3"), mL2.ToString("f3"), Math.Sqrt(mL2).ToString("f2")));
            str.Add(string.Format("{0,20}{1,15}{2,10}{3,10}{4,10}", "S", PVV2.ToString("f3"),
                Sri_sum.ToString("f3"), mS2.ToString("f3"), Math.Sqrt(mS2).ToString("f2")));
            str.Add(string.Format("{0,20}{1,15}", "T", PVVt.ToString("f3")));
            str.Add(string.Format("---------------------------------------------------------------------"));
            if (popt.AdjOpt == ADJOPT.SemiParam)
                str.Add(string.Format("{0,20}{1,-10}", "RegularParam_of_SemiParam:",reg_alpha.ToString("f2")));
            if (popt.WgtDef == Weight.DataSnooping)
            {
                str.Add("\n");
                str.Add(string.Format("{0}{1,-10}", PadLeftEx("数据探测法剔除粗差数:", num), dsp.names.Count));
                str.Add(string.Format("---------------------------------------------------------------------"));
                foreach (DspN item in dsp.names)
                {
                    str.Add(string.Format("{0,20}{1,20}{2,20}", item.type, item.beg, item.end));
                }
                str.Add(string.Format("---------------------------------------------------------------------"));
            }
            if (popt.WgtDef == Weight.Robust || popt.WgtDef == Weight.RobustHelmert)
            {
                str.Add("\n");
                str.Add(string.Format("{0,30}", "抗差估计结果"));
                str.Add(string.Format("---------------------------------------------------------------------"));
                str.Add(string.Format("{0}{1,-10}", PadLeftEx("不变权数:", num), n0));
                str.Add(string.Format("{0}{1,-10}", PadLeftEx("降权数:", num), n1));
                str.Add(string.Format("{0}{1,-10}", PadLeftEx("剔除数:", num), n2));
                str.Add(string.Format("---------------------------------------------------------------------"));
            }
            str.Add("\n");
            File.AppendAllLines(outpath, str);
        }
        //输出方差分量估计结果
        protected void OutputVCE(string outpath)
        {
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30}", "方差分量估计计算"));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0,30}", "单位权方差之比 D(T):D(S)"));
            for (int i = 0; i < cov_res.Count; i++)
            {
                string vce = (cov_res[i].qs / cov_res[i].ql).ToString("f3");
                str.Add(string.Format("{0,20}{1,5}:{2,-10}", cov_res[i].count, 1, vce));
            }
            str.Add(string.Format("---------------------------------------------------------------------"));
            File.AppendAllLines(outpath, str);
        }
        //输出近似坐标 <"append">'true'表示追加，'false'表示新建文件
        protected void OutputApprox(string outpath, bool append)
        {
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30} ", "概略坐标"));
            str.Add(string.Format("----------------------------------------------------------------------"));
            str.Add(string.Format("{0}{1,20}{2,20}",PadLeftEx("Name",sn), "X/m", "Y/m"));
            foreach (KeyValuePair<string, Point2D> kp in knownPoint)
            {
                string X = Approx[kp.Key].X.ToString("f5");
                string Y = Approx[kp.Key].Y.ToString("f5");
                str.Add(string.Format("{0}{1,20}{2,20}",PadLeftEx(kp.Key,sn), X, Y));
            }
            for(int i = 0; i < CopyApprox.Length; i++)
            {
                string name = CopyApprox[i].name;
                if (knownPoint.ContainsKey(name)) continue;
                string x = CopyApprox[i].x.ToString("f5");
                string y = CopyApprox[i].y.ToString("f5");
                str.Add(string.Format("{0}{1,20}{2,20}", PadLeftEx(name, sn), x, y));
            }
            str.Add(string.Format("----------------------------------------------------------------------"));
            str.Add(string.Format("\n"));
            if (append)
                File.AppendAllLines(outpath, str);
            else
                File.WriteAllLines(outpath, str);
        }
        protected double Lri_sum = 0;
        //输出方向平差
        protected void OutputL(string outpath)
        {
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30} ", "方向平差结果"));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0}{1}{2,13}{3,8}{4,8}{5,13}{6,7}",
               PadLeftEx( "From",sn),PadLeftEx("To",sn), "Value/dms", "M/\"", "V/\"", "Result/dms", "Ri"));
            int ri = 0;//第几个方向观测值
            int i = 0;
            var Q1 = P1.Inverse();
            double max = -999, min = 999;
            string MaxSta = "", MaxObs = "", MinSta = "", MinObs = "";
            dsp.names = new List<DspN>();
            foreach (KeyValuePair<string, Dictionary<string, double>> sta in obsL)
            {
                foreach (KeyValuePair<string, double> obs in sta.Value)
                {
                    int ids = FindStaName(sta.Key);
                    int ido = FindObsName(ids, obs.Key);
                    var data = AllObs[ids].data[ido];
                    if (global_popt.WgtDef == Weight.DataSnooping)
                    {
                        if (ri<dsp.idxL.Length && dsp.idxL[ri] != i)
                        {
                            DspN d = new DspN();
                            d.type = "L";
                            d.beg=sta.Key;
                            d.end = obs.Key;
                            dsp.names.Add(d);
                            i++; continue;
                        }
                    }
                    double ML = fai * Math.Sqrt(data.qL);
                    double Lv = D2DMS(DMS2D(obs.Value) + data.vl / 3600);
                    double Ri =data.rocL;
                    str.Add(string.Format("{0}{1}{2,13}{3,8}{4,8}{5,13}{6,7}",
                        PadLeftEx(sta.Key,sn),PadLeftEx(obs.Key,sn), 
                        obs.Value.ToString("f6"), ML.ToString("f2"),
                        data.vl.ToString("f2"), Lv.ToString("f6"), Ri.ToString("f2")));
                    ri++;
                    if (max < Ri)
                    {
                        max = Ri;
                        MaxSta = sta.Key;
                        MaxObs = obs.Key;
                    }
                    if (min > Ri)
                    {
                        min = Ri;
                        MinSta = sta.Key;
                        MinObs = obs.Key;
                    }
                    Lri_sum += Ri;
                    i++;
                }
            }
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("{0,30}{1,-20}", "方向最小多余观测分量: ", min.ToString("f2") + " (" + MinSta + "-->" + MinObs + ")"));
            str.Add(string.Format("{0,30}{1,-20}", "方向最大多余观测分量: ", max.ToString("f2") + " (" + MaxSta + "-->" + MaxObs + ")"));
            str.Add(string.Format("{0,30}{1,-20}", "方向平均多余观测分量: ", (Lri_sum / nl).ToString("f2")));
            str.Add(string.Format("{0,30}{1,-20}", "方向多余分量总和: ", Lri_sum.ToString("f2")));
            str.Add(string.Format("---------------------------------------------------------------------"));
            str.Add(string.Format("\n"));
            File.AppendAllLines(outpath, str);
        }
        protected double Sri_sum = 0;
        //输出距离平差
        protected void OutputS(string outpath)
        {
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30}", "距离平差结果"));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0}{1}{2,13}{3,8}{4,8}{5,13}{6,7}",
               PadLeftEx( "From",sn),PadLeftEx( "To",sn), "Value/m", "M/mm", "V/mm", "Result/m", "Ri"));
            int ri = 0;//第几个距离观测值
            int t = V1.RowCount;//L个数
            int i = 0;
            var Q2 = P2.Inverse();
            double max = -999, min = 999;
            string MaxSta = "", MaxObs = "", MinSta = "", MinObs = "";
            foreach (KeyValuePair<string, Dictionary<string, double>> sta in obsS)
            {
                foreach (KeyValuePair<string, double> obs in sta.Value)
                {
                    int ids = FindStaName(sta.Key);
                    int ido = FindObsName(ids, obs.Key);
                    var data = AllObs[ids].data[ido];
                    if (global_popt.WgtDef == Weight.DataSnooping)
                    {
                        if (dsp.idxS.Length != 0 && dsp.idxS[ri] != i)
                        {
                            DspN d = new DspN();
                            d.type = "S";
                            d.beg = sta.Key;
                            d.end = obs.Key;
                            dsp.names.Add(d);
                            i++; continue;
                        }
                    }
                    double MS = fai * Math.Sqrt(data.qS);
                    double Sv = obs.Value + data.vs / 1000;
                    double Ri = data.rocS;
                    str.Add(string.Format("{0}{1}{2,13}{3,8}{4,8}{5,13}{6,7}",
                        PadLeftEx( sta.Key,sn), PadLeftEx(obs.Key, sn), 
                        obs.Value.ToString("f5"), MS.ToString("f2"),
                        data.vs.ToString("f2"), Sv.ToString("f5"), Ri.ToString("f2")));
                    ri++;
                    if (max < Ri)
                    {
                        max = Ri;
                        MaxSta = sta.Key;
                        MaxObs = obs.Key;
                    }
                    if (min > Ri)
                    {
                        min = Ri;
                        MinSta = sta.Key;
                        MinObs = obs.Key;
                    }
                    Sri_sum += Ri;
                    i++;
                }
            }
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0,30}{1,-20}", "边长最小多余观测分量: ", min.ToString("f2") + " (" + MinSta + "-->" + MinObs + ")"));
            str.Add(string.Format("{0,30}{1,-20}", "边长最大多余观测分量: ", max.ToString("f2") + " (" + MaxSta + "-->" + MaxObs + ")"));
            str.Add(string.Format("{0,30}{1,-20}", "边长平均多余观测分量: ", (Sri_sum / ns).ToString("f2")));
            str.Add(string.Format("{0,30}{1,-20}", "边长多余分量总和: ", Sri_sum.ToString("f2")));
            str.Add("---------------------------------------------------------------------");
            str.Add("\n");
            File.AppendAllLines(outpath, str);
        }
        //输出点位信息
        protected void OutputXY(string outpath, pindex pIndex)
        {
            if (Result.Count < PointName.Count) return;
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30}", "平差坐标及其精度"));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0}{1,16}{2,15}{3,8}{4,8}{5,8}{6,8}{7,8}{8,11}",
               PadLeftEx("Name",sn), "X/m", "Y/m", "MX/mm", "MY/mm", "MP/mm", "E/mm", "F/mm", "T/dms"));
            foreach (KeyValuePair<string, Point2D> p in knownPoint)
            { //先输出已知点
                string name = p.Key;
                string X = p.Value.X.ToString("f5");
                string Y = p.Value.Y.ToString("f5");
                str.Add(string.Format("{0}{1,16}{2,15}",PadLeftEx( p.Key,13), X, Y));
            }
            string weakP = string.Empty;
            double weakX = 0, weakY = 0, weakMX = 0, weakMY = 0, weakMP = 0,
                weakE = 0, weakF = 0, weakT = 0;
            double MxSum = 0, MySum = 0, MpSum = 0;
            foreach (KeyValuePair<string, int> p in pIndex)
            {//未知点
                double Te = 0, T = 0;
                double X = Result[p.Key].X;
                double Y = Result[p.Key].Y;
                //*****************计算坐标协因数******************
                double Qxx = QX[2 * p.Value, 2 * p.Value];
                double Qyy = QX[2 * p.Value + 1, 2 * p.Value + 1];
                double Qxy = QX[2 * p.Value, 2 * p.Value + 1];
                double MX = sigma * Math.Sqrt(Qxx);
                double MY = sigma * Math.Sqrt(Qyy);
                double MP = Math.Sqrt(MX * MX + MY * MY);
                MxSum += MX; MySum += MY; MpSum += MP;
                //*****************计算误差椭圆******************
                double K = Math.Sqrt(Math.Pow(Qxx - Qyy, 2) + 4 * Qxy * Qxy);
                double E = sigma * Math.Sqrt(0.5 * (Qxx + Qyy + K));
                double F = sigma * Math.Sqrt(0.5 * (Qxx + Qyy - K));
                //*************误差椭圆极大方向①*******
                double Qe = 0.5 * (Qxx + Qyy + K);
                Te = Keep_0To360(RD2D(Math.Atan((Qe - Qxx) / Qxy)));
                if (Te >= 180) //极大方向是反方向的两个，取较小的
                    Te -= 180;
                T = D2DMS(Te);
                Result[p.Key].Mx = MX;
                Result[p.Key].My = MY;
                Result[p.Key].Mp = MP;
                Result[p.Key].E = E;
                Result[p.Key].F = F;
                Result[p.Key].T = T;
                ////**********②***********
                //double TE = CorrectD(RD2D(0.5 * Math.Atan(2 * Qxy / (Qxx - Qyy))));
                //if (Qxy > 0)
                //{
                //    if (!Location13(TE))
                //        TE += 90;
                //}
                //else
                //{
                //    if (Location13(TE))
                //        TE += 90;
                //}
                //TE = CorrectD(TE);
                //if (TE >= 180) //极大方向是反方向的两个，取较小的
                //    TE -= 180;
                //Error0.Add(TE);
                //double T = D2DMS(TE);
                ////***********************************

                if (weakMP < MP)
                { //寻找最弱点
                    weakMP = MP; weakMX = MX; weakMY = MY; weakP = p.Key;
                    weakX = X; weakY = Y; weakE = E; weakF = F; weakT = T;
                }
                str.Add(string.Format("{0}{1,16}{2,15}{3,8}{4,8}{5,8}{6,8}{7,8}{8,11}",PadLeftEx(p.Key,sn),
                  X.ToString("f5"), Y.ToString("f5"), MX.ToString("f2"), MY.ToString("f2"),
                  MP.ToString("f2"), E.ToString("f2"), F.ToString("f2"), T.ToString("f4")));
            }
            str.Add("---------------------------------------------------------------------");
            string mx_ave = (MxSum / ukp_num).ToString("f2");
            string my_ave = (MySum / ukp_num).ToString("f2");
            string mp_ave = (MpSum / ukp_num).ToString("f2");
            str.Add(string.Format("{0,15}{1,-10}{2,15}{3,-10}{4,15}{5,-10}",
                "MX均值(mm):", mx_ave, "MY均值(mm):", my_ave, "MP均值(mm):", mp_ave));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0,30}", "最弱点及其精度"));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0}{1,16}{2,15}{3,8}{4,8}{5,8}{6,8}{7,8}{8,11}",
               PadLeftEx( "Name",sn), "X/m", "Y/m", "MX/mm", "MY/mm", "MP/mm", "E/mm", "F/mm", "T/dms"));
            str.Add(string.Format("{0}{1,16}{2,15}{3,8}{4,8}{5,8}{6,8}{7,8}{8,11}",PadLeftEx(weakP,sn),
                weakX.ToString("f5"), weakY.ToString("f5"), weakMX.ToString("f2"), weakMY.ToString("f2"),
                weakMP.ToString("f2"), weakE.ToString("f2"), weakF.ToString("f2"), weakT.ToString("f4")));
            str.Add("---------------------------------------------------------------------");
            str.Add("\n");
            File.AppendAllLines(outpath, str);

#if false    // 输出陀螺自适应定权对最弱点的影响
            Console.WriteLine(string.Format("{0,5},{1,5},{2,5},{3,5}", "Adp="+global_popt.GyOuAdp, 
                "Name="+weakP,"横坐标="+weakX.ToString("f3"),"M="+weakE.ToString("f0")));
#endif
        }
        //输出相对点位
        protected void OutputLS(string outpath, pindex pIndex)
        {
            double sumMA = 0, sumMS = 0;
            List<string> str = new List<string>();
            str.Add(string.Format("{0,30}", "方位角、边长及其相对精度"));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0}{1}{2,13}{3,8}{4,13}{5,8}{6,9}{7,8}{8,8}{9,11}",
               PadLeftEx("From", sn), PadLeftEx("To", sn), "Azi/dms", "MA/\"", "S/m", "MS/mm", "S/MS", "E/mm", "F/mm", "T/dms"));
            //var QX = N.Inverse();//坐标协因数
            int count = 0; //第几对相对点位，第几个照准点
            int ri = 0;
            double weakAzi = 0;
            double weakS_MS = 0; string weakSMS = "";//最弱边
            string weakLi = "", weakLj = "", weakMA = "";
            string weakS = "", weakMS = "";
            string weakEE = "", weakFF = "", weakTT = "";
            foreach (KeyValuePair<string, Dictionary<string, double>> sta in obsL)
            {
                foreach (KeyValuePair<string, double> obs in sta.Value)
                {
                    if (global_popt.WgtDef== Weight.DataSnooping)
                    {
                        if (dsp.idxL.Length != 0 && dsp.idxL[ri] != count)
                        {
                            count++; continue;
                        }
                    }
                    double dx = Result[obs.Key].X - Result[sta.Key].X;
                    double dy = Result[obs.Key].Y - Result[sta.Key].Y;
                    double S2 = dx * dx + dy * dy;
                    double azi = RD2DMS(GetAzi(dx, dy));
                    //*************区分已知点和未知点************
                    int i = 0, j = 0, I = 0, J = 0;
                    if (pIndex.ContainsKey(sta.Key))
                    {
                        i = pIndex[sta.Key];
                        I = 1;
                    }
                    else I = 0;
                    if (pIndex.ContainsKey(obs.Key))
                    {
                        j = pIndex[obs.Key];
                        J = 1;
                    }
                    else J = 0;  //若j点是已知点，Qdx=Qxi,Qdy=Qyi,Qdxy=Qxi.
                    if (I == 0 && J == 0)
                        continue;
                    var Qxixi = QX[2 * i, 2 * i]; var Qxjxj = QX[2 * j, 2 * j];
                    var Qxixj = QX[2 * i, 2 * j]; var Qxjxi = Qxixj;
                    var Qyiyi = QX[2 * i + 1, 2 * i + 1]; var Qyjyj = QX[2 * j + 1, 2 * j + 1];
                    var Qyiyj = QX[2 * i + 1, 2 * j + 1]; var Qyjyi = Qyiyj;
                    var Qxiyi = QX[2 * i, 2 * i + 1]; var Qyixi = Qxiyi;
                    var Qxjyj = QX[2 * j, 2 * j + 1]; var Qyjxj = Qxjyj;
                    var Qxiyj = QX[2 * i, 2 * j + 1]; var Qyjxi = Qxiyj;
                    var Qxjyi = QX[2 * j, 2 * i + 1]; var Qyixj = Qxjyi;
                    //***************方位角、距离及其精度************
                    double[,] qxy = {{ I*Qxixi,     I*Qxiyi,     I*J*Qxixj, I*J*Qxiyj},
                                     { I*Qyixi,     I*Qyiyi,     I*J*Qyixj, I*J*Qyiyj},
                                     { I*J*Qxjxi,   I*J*Qxjyi,   J*Qxjxj ,  J*Qxjyj  },
                                     { I*J*Qyjxi,   I*J*Qyjyi,   J*Qyjxj,   J*Qyjyj  } };
                    double a = rou * dy / S2 / 1000;
                    double b = rou * dx / S2 / 1000;
                    double[,] ka = { { a, -1 * b, -1 * a, b } };
                    var QXY = mat.Build.DenseOfArray(qxy);
                    var KA = mat.Build.DenseOfArray(ka);
                    var QA = KA * QXY * KA.Transpose();
                    double MA = sigma * Math.Sqrt(QA[0, 0]);
                    sumMA += MA;
                    double S = Math.Sqrt(dx * dx + dy * dy);
                    double c = dx / S; double d = dy / S;
                    double[,] ks = { { -1 * c, -1 * d, c, d } };
                    var KS = mat.Build.DenseOfArray(ks);
                    var QS = KS * QXY * KS.Transpose();
                    double MS = sigma * Math.Sqrt(Math.Abs(QS[0, 0]));
                    sumMS += MS;
                    double sms = S * 1000 / MS;
                    int S_MS = Convert.ToInt32(sms);
                    //*********计算相对误差椭圆********//
                    double Qdx = I * Qxixi + J * Qxjxj - I * J * 2 * Qxixj;
                    double Qdy = I * Qyiyi + J * Qyjyj - I * J * 2 * Qyiyj;
                    double Qdxy = I * Qxiyi - I * J * Qxiyj - I * J * Qxjyi + J * Qxjyj;
                    double K = Math.Sqrt(Math.Pow(Qdx - Qdy, 2) + 4 * Qdxy * Qdxy);
                    double E = sigma * Math.Sqrt(0.5 * (Qdx + Qdy + K));
                    double F = sigma * Math.Sqrt(0.5 * (Qdx + Qdy - K));
                    double Qe = 0.5 * (Qdx + Qdy + K);
                    double Te = Keep_0To360(RD2D(Math.Atan((Qe - Qdx) / Qdxy)));
                    if (Te >= 180) //极大方向是反方向的两个，取较小的
                        Te -= 180;
                    string T = D2DMS(Te).ToString("f4");
                    string SMS = S_MS.ToString();
                    //*********************计算最弱边********************
                    str.Add(string.Format("{0}{1}{2,13}{3,8}{4,13}{5,8}{6,9}{7,8}{8,8}{9,11}",
                       PadLeftEx( sta.Key,sn),PadLeftEx( obs.Key,sn), 
                       azi.ToString("f6"), MA.ToString("f2"), 
                       S.ToString("f5"),MS.ToString("f2"), SMS, 
                        E.ToString("f2"), F.ToString("f2"), T));
                    if (count == 0)
                    { //最弱边赋初值
                        weakLi = sta.Key; weakLj = obs.Key;
                        weakAzi = azi; weakMA = MA.ToString("f2");
                        weakS = S.ToString("f5"); weakMS = MS.ToString("f2"); weakS_MS = S_MS; weakSMS = SMS;
                        weakEE = E.ToString("f2"); weakFF = F.ToString("f2"); weakTT = T;
                    }
                    else
                    {
                        if (weakS_MS > S_MS)
                        {//边长相对中误差的倒数越小，边长相对中误差就越大，精度越弱
                            weakLi = sta.Key; weakLj = obs.Key;
                            weakAzi = azi; weakMA = MA.ToString("f2");
                            weakS = S.ToString("f5"); weakMS = MS.ToString("f2"); weakS_MS = S_MS; weakSMS = SMS;
                            weakEE = E.ToString("f2"); weakFF = F.ToString("f2"); weakTT = T;
                        }
                    }
                    count++;ri++;
                }
            }
            str.Add(string.Format("----------------------------------------------------------------------"));
            string aveMA = (sumMA / nl).ToString("f2");
            string aveMS = (sumMS / nl).ToString("f2");
            str.Add(string.Format("{0,20}{1,-10}{2,10}{3,-10}", "MA均值:", aveMA, "MS均值:", aveMS));
            str.Add(string.Format("----------------------------------------------------------------------"));
            str.Add(string.Format("{0,30}", "最弱边及其相对精度"));
            str.Add("---------------------------------------------------------------------");
            str.Add(string.Format("{0}{1}{2,13}{3,8}{4,13}{5,8}{6,9}{7,8}{8,8}{9,11}",
               PadLeftEx("From", sn), PadLeftEx("To", sn), "Azi/dms", "MA/\"", "S/m", "MS/mm", "S/MS", "E/mm", "F/mm", "T/dms"));
            str.Add(string.Format("{0}{1}{2,13}{3,8}{4,13}{5,8}{6,9}{7,8}{8,8}{9,11}",PadLeftEx( weakLi,sn),
            PadLeftEx( weakLj,sn), weakAzi.ToString("f6"), weakMA, weakS, weakMS, weakSMS, weakEE, weakFF, weakTT));
            str.Add("---------------------------------------------------------------------");
            str.Add("\n");
            File.AppendAllLines(outpath, str);
        }

        //输出ou2文件
        protected void Output(Proc_t popt, string outpath)
        {
            OutputTitle(popt, outpath);//标题
            OutputApprox(outpath, true);//概略坐标
            if (nl != 0)
                OutputL(outpath); //L
            if (ns != 0)
                OutputS(outpath);//S
            OutputXY(outpath, ukpIndex);//点位
            OutputLS(outpath, ukpIndex);//相对点位
            if (popt.GyOpt.UseGyro == 1)
                OutputGyro(outpath);
            OutputProfile(popt, outpath);//概况
            if (popt.WgtDef >= Weight.Helmert && popt.WgtDef < Weight.Robust || popt.WgtDef == Weight.RobustHelmert)
                OutputVCE(outpath);//定权
        }

        // 绘图文件输出
        void Plot(string path)
        {
            List<string> str = new List<string>();
            str.Add(string.Format("{0,10}{1,10}", Approx.Count, n));
            foreach (KeyValuePair<string, Point2D> kp in knownPoint)
            {
                string X = Approx[kp.Key].X.ToString("f5");
                string Y = Approx[kp.Key].Y.ToString("f5");
                str.Add(string.Format("{0}{1,16}{2,16}", PadLeftEx(kp.Key, sn), X, Y));
            }
            foreach (KeyValuePair<string, int> p in ukpIndex)
            {
                string X = Result[p.Key].X.ToString("f5");
                string Y = Result[p.Key].Y.ToString("f5");
                string E = (Result[p.Key].E * 0.1).ToString("f2");
                string F = (Result[p.Key].F * 0.1).ToString("f2");
                string T = Result[p.Key].T.ToString("f4");
                str.Add(string.Format("{0}{1,16}{2,16}{3,8}{4,8}{5,13}",
                    PadLeftEx(p.Key, sn), X, Y, E, F, T));
            }
            foreach (KeyValuePair<string, Dictionary<string, double>> sta in obsL)
            {
                foreach (KeyValuePair<string, double> obs in sta.Value)
                {
                    str.Add(string.Format("{0}{1}{2,10}",
                        PadLeftEx(sta.Key, sn), PadLeftEx(obs.Key, sn), "L"));
                }
            }
            foreach (KeyValuePair<string, Dictionary<string, double>> sta in obsS)
            {
                foreach (KeyValuePair<string, double> obs in sta.Value)
                {
                    str.Add(string.Format("{0}{1}{2,10}",
                        PadLeftEx(sta.Key, sn), PadLeftEx(obs.Key, sn), "S"));
                }
            }
            File.WriteAllLines(path, str);
        }
    }
}
